1. Stats de base, modèles linéaires, GLM

Distribution statistique

Décrit les probabilités d’apparition des valeurs d’une distribution donnée.

Distribution normale

On parle aussi souvent d’une distribution Gaussienne.


\[ f(x) = \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2} \]


\(\mu\) et \(\sigma\) sont des paramètres, respectivement la moyenne et l’écart-type.

Forme de la distribution normale

x<-seq(-5,15,by=0.01)
mu<-5
sigma<-2
y<-(1/(sigma*sqrt(2*pi))*exp((-1/2)*((x-mu)/sigma)^2))
plot(x,y,type="l")

Simuler une distribution normale

On peut utiliser la fonction rnorm pour simuler des valeurs provenant d’une distribution normale.

set.seed(123)
n<-1000
x<-rnorm(n,mu,sigma)
hist(x,breaks=50)

Interprétation

v<-seq(-5,15,by=0.01)
hist(x,breaks=50,freq=FALSE,xlim=range(v))
curve(dnorm(x,mu,sigma),add=TRUE)
y<-dnorm(v,mu,sigma)
polygon(c(v[v>=6 & v<=8],8,6), c(y[v>=6 & v<=8],0,0),col=alpha("red",0.5),border=NA)

L’aire sous la courbe est de 1 et la proportion de la surface en rouge représente la probabilité d’obtenir une valeur entre 6 et 8 avec une moyenne de \(\mu = 5\) et un écart-type de \(\sigma = 2\). On peut estimer cette probabilité avec la fonction pnorm intégrée à R. Ce type de fonction est intégré pour plusieurs distribution statistique (voir ?Distributions).

Probabilité d’obtenir une valeur?
pnorm(5,mean=5,sd=2) # probabilité d'obtenir une valeur <= 5
## [1] 0.5
pnorm(-2,mean=5,sd=2) # probabilité d'obtenir une valeur <= -2
## [1] 0.0002326291
pnorm(4,mean=5,sd=2) # probabilité d'obtenir une valeur <= 4
## [1] 0.3085375

Probabilité d’obtenir une valeur entre 6 et 8

pnorm(8,5,2)-pnorm(6,5,2)
## [1] 0.2417303

Distribution d’une variable


\[ X \sim \mathcal{N}(\mu,\,\sigma^{2})\,\]


\(X\) suit une distribution normale avec une moyenne \(\mu\) et une variance \(\sigma^{2}\)

Erreur-type


\[ se = \frac{sd}{\sqrt{n}} \]


ou


\[ \sigma_{\bar{x}} = \frac{\sigma_{x}}{\sqrt{n}} \]


n<-100000
pop<-rnorm(n,100,20) # population fictive de n individus avec une moyenne de 100 et écart-type de 20
ech<-sample(pop,30) # échantillon aléatoire de 30 individus
mean(ech) # moyenne
## [1] 100.5529
sd(ech)/sqrt(30) # erreur-type
## [1] 3.27161

Simulation

Ici, on refait cet échantillonnage nreps fois et on calcule l’écart-type des différentes moyennes obtenues.

nreps<-1000
ech1000<-sapply(1:nreps,function(i){
  s<-sample(pop,30) 
  mean(s)
})
hist(ech1000)

Écart-type des valeurs simulées

sd(ech1000)
## [1] 3.734432

Comparons ceci à l’erreur-type calculé à partir de notre seul échantillon:

sd(ech)/sqrt(30) # erreur-type
## [1] 3.27161
Interprétation de l’erreur-type

L’écart-type de ces moyennes correspond à l’erreur-type estimée à partir d’un seul échantillon. En d’autres mots, l’erreur-type représente:

  • l’écart-type si des moyennes obtenues si on refaisait plusieurs fois le même échantillonnage

  • une mesure de précision dans l’estimation de la moyenne


Le concept d’erreur-type s’applique aussi à d’autres paramètres estimés et non seulement à la moyenne.


L’erreur-type d’une statistique ou d’un paramètre est une estimation de l’écart-type de sa distribution d’échantillonnage.

Intervalle de confiance

Reprenons notre population fictive pop qui a une moyenne de 100 et un écart-type de 20. On se rappelle qu’un intervalle de confiance pour une moyenne provenant d’un échantillon distribué normalement (ou presque) se calcule de la façon suivante:


\[\bar{x} ± t_{dl,\alpha/2}\sigma_{\bar{x}}\]


On peut calculer cet intervalle de confiance avec R:

mean(ech)+qt(0.025,df=length(ech)-1,lower.tail=FALSE)*(sd(ech)/sqrt(length(ech)))*c(-1,1)
## [1]  93.86168 107.24407

Simulation

Maintenant, reprenons note population fictive et calculons à chaque fois cet intervalle de confiance. Ensuite, calculons quelle est la proportion des intervalles qui contiennent la véritable moyenne de 100.

nreps<-1000
ech1000<-lapply(1:nreps,function(i){
  s<-sample(pop,30) 
  c(mean(s),mean(s)+qt(0.025,df=length(s)-1,lower.tail=FALSE)*(sd(s)/sqrt(length(s)))*c(-1,1))
})
ci<-as.data.frame(do.call("rbind",ech1000))
names(ci)<-c("mean","lowerCI","upperCI")

sum(100>=ci$lowerCI & 100<=ci$upperCI)/nreps
## [1] 0.952

Interprétation de l’intervalle de confiance

Un intervalle de confiance à 95% veut dire que 95% des intervalles de confiance qu’on obtiendrait si on refaisait le même échantillonnage contiendront la véritable moyenne de la population.

Plus généralement, on devrait parler d’un paramètre (moyenne, variance, pente, etc.) plutôt que d’une moyenne.

On entend souvent à tort que cela veut dire qu’il y a 95% de chances que le paramètre soit à l’intérieur des bornes. Or le paramètre est considéré comme étant fixe et la probabilité porte sur le fait que l’intervalle contienne la véritable valeur du paramètre ou pas.

Tout comme l’erreur-type, l’intervalle de confiance est une mesure de précision dans l’estimation d’un paramètre.

Valeur de p

Un exemple avec le test de \(t\) et une population fictive.

set.seed(1)
n<-10000
pop1<-rnorm(n,20,3) # population fictive de n individus avec une moyenne de 20
pop2<-rnorm(n,21,3) # population fictive de n individus avec une moyenne de 21

ech1<-sample(pop1,10) # échantillon de 10 individus de chaque population
ech2<-sample(pop2,10)

test<-t.test(ech1,ech2)
test
## 
##  Welch Two Sample t-test
## 
## data:  ech1 and ech2
## t = -1.189, df = 14.686, p-value = 0.2533
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -4.030060  1.147318
## sample estimates:
## mean of x mean of y 
##  19.22124  20.66261

Interprétation de la Valeur de p

plot(0,0,xlim=c(-4,4),ylim=c(0,0.5),type="n")
curve(dt(x,test$parameter),add=TRUE)
v<-seq(-5,15,by=0.01)
y<-dt(v,test$parameter)
polygon(c(v[v<=test$statistic],test$statistic,min(v)), c(y[v<=test$statistic],0,0),col=alpha("red",0.5),border=NA)
polygon(c(v[v>=abs(test$statistic)],max(v),abs(test$statistic)), c(y[v>=abs(test$statistic)],0,0),col=alpha("red",0.5),border=NA)
abline(v=rep(test$statistic,2)*c(1,-1),lwd=2,lty=3,col="red")

Simulations

Reprenons le test de \(t\) précédent, mais cette fois, prenons nos échantillons dans la même population (pop1) à chaque fois.

nreps<-1000

pvalue<-sapply(1:nreps,function(i){
  ech1<-sample(pop1,5) 
  ech2<-sample(pop1,5)
  test<-t.test(ech1,ech2)
  test$p.value#statistic
})

hist(pvalue,breaks=50)
abline(v=0.05,col="red",lwd=2,lty=3)

sum(pvalue<=0.05)/length(pvalue)
## [1] 0.051

Modèles linéaires simples


On a ici la formulation classique ou habituelle d’un modèle linéaire simple.


\[ y = \beta_{0} + \beta_{1}x_{1} + \beta_{2}x_{2} + ... + \beta_{n}x_{n} + \epsilon \]

\[ \epsilon \sim \mathcal{N}(0,\,\sigma^{2})\ \]

Autre formulation


\[ \mu = \beta_{0} + \beta_{1}x_{1} + \beta_{2}x_{2} + ... + \beta_{n}x_{n}\]

\[ y \sim \mathcal{N}(\mu,\,\sigma^{2})\ \]



\(y\) = observations

\(\mu\) = prédicteur linéaire

On a donc un modèle représentant le lien entre nos variables et on cherche à estimer les paramètres de ce modèles à l’aide de données.

Suppositions de bases

Quelles sont les suppositions de bases d’un modèle linéaire tel que celui présenté?

  • Normalité des résidus?

  • Homoscédasticité de la variance (la variance ne dépend pas de ou ne varie pas avec \(x_{n}\))

  • Linéarité des relations entre \(x_{n}\) et \(y\)?

  • Corrélation entre les \(x_{n}\)?

  • Indépendence entre les observations?

Simulations

  • Une des meilleures façons de s’assurer qu’on comprend comment fonctionne un modèle

  • Permet d’étudier le comportement d’un modèle avec des valeurs connues

Simulons un modèle

# nb d'observations
n<-75

# valeurs fictives des variables explicatives selon une distribution uniforme
x1<-runif(n,0,100)
x2<-runif(n,0,10)
x3<-runif(n,0,1)

# valeurs des paramètres beta
b0<-100
b1<-0.2
b2<-1
b3<--10

# prédicteur linéaire
mu<-b0+b1*x1+b2*x2+b3*x3

# observations à partir d'une distribution normale avec un écart-type de 5
y<-rnorm(n,mu,5)

# données
d<-data.frame(y,x1,x2,x3)

Fonction lm

La fonction lm permet de faire des modèles linéaires simples, ce qui inclue les régressions simples, les ANOVAs, ANCOVAs, régression multiples, etc.

Remarquez l’utilisation de l’argument data où la fonction ira chercher les différents éléments spécifiés dans la formule.

m<-lm(y~x1+x2+x3,data=d)

Coefficients

On peut appliquer plusieurs fonctions pour extraire les éléments d’un modèle, en particulier un modèle effectué avec la fonction lm.

coef(m)
## (Intercept)          x1          x2          x3 
##  96.9666068   0.1882304   1.0717263  -4.0315717

Sommaire

La fonction summary est la plus utile pour extraire l’information d’un modèle.

summary(m)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.518  -2.927   0.206   2.739  12.382 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 96.96661    1.80048  53.856  < 2e-16 ***
## x1           0.18823    0.01909   9.858 6.23e-15 ***
## x2           1.07173    0.18827   5.693 2.62e-07 ***
## x3          -4.03157    1.92616  -2.093   0.0399 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.02 on 71 degrees of freedom
## Multiple R-squared:  0.6444, Adjusted R-squared:  0.6294 
## F-statistic: 42.89 on 3 and 71 DF,  p-value: 6.301e-16

Vérification des suppositions de bases

On peut utiliser la fonction plot qui lorsqu’elle est appliquée sur un objet de classe lm produira des graphiques permettant de vérifier certaines suppositions de bases ou conditions d’application.

par(mfrow=c(2,2))
plot(m)

Fonctions accessoires

On peut également créer facilement certains de ces graphiques avec d’autres fonctions.

par(mfrow=c(1,3))
plot(resid(m),fitted(m))
hist(resid(m))
qqnorm(resid(m))
qqline(resid(m))

Facteurs

# simulations d'un facteur
d$fac<-factor(sample(c("A","B","C"),nrow(d),replace=TRUE))
d$y<-d$y+2*as.integer(d$fac)-1

# modèle modifié avec le facteur
m<-lm(y~x1+x2+x3+fac,data=d)
summary(m)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + fac, data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.3549  -3.0036   0.3166   2.6305  13.2288 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  97.2099     1.9870  48.924  < 2e-16 ***
## x1            0.1861     0.0190   9.794 1.09e-14 ***
## x2            1.0464     0.1854   5.645 3.39e-07 ***
## x3           -3.8057     1.8956  -2.008   0.0486 *  
## facB          1.9366     1.4416   1.343   0.1836    
## facC          6.5096     1.3744   4.736 1.12e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.932 on 69 degrees of freedom
## Multiple R-squared:  0.696,  Adjusted R-squared:  0.674 
## F-statistic:  31.6 on 5 and 69 DF,  p-value: < 2.2e-16

Comprendre son modèle avec des graphiques

La fonction predict sert à calculer la valeur prédite par le modèle pour chaque observation. Par exemple, ceci permet de mettre en relation les valeurs observées et les valeurs prédites. La fonction predict dans ce cas calcule la valeur prédite par le modèle pour chaque observation dans la base de donnée.

plot(d$y,predict(m))

Soumettre de nouvelle valeurs à predict

Plus souvent, on veut illustrer l’effet de nos variables et pour cela, il faut soumettre des valeurs à la fonction predict. Cette fonction est très utile lorsque l’on veut un maximum de contrôle sur les valeurs à prédire.

x<-seq(min(d$x1),max(d$x1),length.out=50)
nd<-data.frame(x1=x,x2=mean(d$x2),x3=mean(d$x3),fac="B") # on fixe les autres variables à leur valeur moyenne
p<-predict(m,newdata=nd)
plot(d$x1,d$y) # observations
lines(x,p) # valeurs prédites

Package visreg

Le package visreg permet d’illustrer facilement les prédictions d’un modèle (ou les effets marginaux des différentes variables explicatives). Il utilise en arrière-plan la fonction predict.

library(visreg)
par(mfrow=c(2,2))
visreg(m)

Package ggeffects

Le package ggeffects permet également d’illustrer les prédictions d’un modèle, mais ce package se base sur l’utilisation du package ggplot2 pour construire les graphiques. Le principe demeure le même: en arrière-plan la fonction predict est utilisée. On peut également combiner les différents graphiques générés en utilisant le package patchwork et sa fonction wrap_plots.

library(ggeffects)
library(patchwork)
g<-ggpredict(m)
p<-plot(g,add.data=TRUE)
wrap_plots(p)

Interactions

Dans vos mots, qu’est-ce qu’une interaction?


C’est lorsque l’effet d’une variable dépend des valeurs d’une autre variable (ou de plusieurs).

Un exemple

# modèle modifié avec le facteur
m<-lm(y~x1+x2+x3+x1*fac,data=d)
summary(m)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x1 * fac, data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.0406  -1.9485   0.1316   2.1335  14.2374 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 98.50371    2.15978  45.608  < 2e-16 ***
## x1           0.14204    0.03151   4.508  2.7e-05 ***
## x2           1.15443    0.18388   6.278  2.9e-08 ***
## x3          -2.44076    1.86546  -1.308  0.19521    
## facB         0.66407    2.63788   0.252  0.80201    
## facC        -1.01455    3.01397  -0.337  0.73746    
## x1:facB      0.01671    0.04332   0.386  0.70090    
## x1:facC      0.13612    0.04871   2.794  0.00678 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.692 on 67 degrees of freedom
## Multiple R-squared:  0.7329, Adjusted R-squared:  0.705 
## F-statistic: 26.26 on 7 and 67 DF,  p-value: < 2.2e-16

Illustration graphique

visreg(m,"x1",by="fac",overlay=TRUE)


On a une interaction significative, mais est-elle importante d’un point de vue biologique? Peut-être pas… Toujours distinguer entre la signification statistique et la signification biologique.

Une valeur de p significative c’est bien beau, mais ce qui importe vraiment c’est la taille de l’effet.

Interactions avec ggeffects
g<-ggpredict(m,terms=c("x1","fac"),message=FALSE)
plot(g,add.data=TRUE,jitter=FALSE)

Interprétion

  • Beaucoup plus facile avec un graphique qu’en regardant une table de coefficients.

  • En présence d’une interaction (significative), les effets simples sont plus difficilement interprétables et ne peuvent être interprétés sans faire référence à l’interaction (à moins d’ajustements particuliers voir Schielzeth 2010).

  • Si pour deux variables impliquées dans une interaction les effets simples ne sont pas significatifs, mais que leur interaction l’est, il faut considérer que ces deux variables ont un effet même si les coefficients associés aux effets simples ne sont pas significatifs.

Types d’erreurs


  • Erreur de type I: rejeter l’hypothèse nulle, alors qu’elle est vraie (faux positif).


  • Erreur de type II: ne pas rejeter l’hypothèse nulle, alors qu’elle est fausse (faux négatif).

Simulations

Simulons un modèle dans lequel les variables n’ont aucun effet.

Que se passera-t-il si on répète plusieurs fois ce scénario?

À quoi ressembleront nos valeurs de p?

Créons d’abord une fonction sim_model pour générer un modèle linéaire fictif.

sim_model<-function(b0,b1,b2,b3){
  n<-200
  x1<-runif(n,0,100)
  x2<-runif(n,0,10)
  x3<-runif(n,0,1)
  mu<-b0+b1*x1+b2*x2+b3*x3
  y<-rnorm(n,mu,5)
  d<-data.frame(y,x1,x2,x3)
  m<-lm(y~x1+x2+x3,data=d)
  m
}

Simulons

Ensuite, simulons nsims modèles fictifs et emmagasinons ces modèles dans une liste nommée list_models.

nsims<-500

list_models<-lapply(1:nsims,function(i){
  sim_model(b0=100,b1=0,b2=0,b3=0)  
})


Si on illustre les valeurs de p obtenues pour chacun des coefficients associés aux variables \(x_{n}\) à l’aide d’un histogramme, à quoi ressemblera cet histogramme?

Valeurs de p obtenus

p<-lapply(list_models,function(i){
  summary(i)$coef[2:4,4]
})

hist(unlist(p),breaks=seq(0,1,by=0.05),xlab="Valeurs de p")

Conclusions
  • Dans une certaine proportion des cas ( ~ 5% ), on conclut à un effet, alors qu’il n’y en a pas pas! C’est le fameux seuil \(\alpha\) de 0.05.

  • Sous l’hypothèse nulle, les valeurs de p sont distribuées uniformément entre 0 et 1.

Modèles linéaires standards

\[y \sim \mathcal{N}(\mu, \sigma^2)\]

\[\mu=\beta_0 + \beta_1 X_1 + \beta_2 X_2 + ... + \beta_k X_k\]

Formulation matricielle

\[y \sim \mathcal{N}(\beta\mathbf{X}, \mathbf{I}\sigma^2)\]


\(\mathbf{X}\) est une matrice avec les valeurs des variables explicatives.

\(\mathbf{I}\) est une matrice identité qui a pour effet d’assigner une variance égale \(\sigma^2\) à toutes les observations.

Modèles linéaires généraux

\[y \sim \mathcal{N}(\beta\mathbf{X}, \mathbf{\Sigma})\]

\(\mathbf{\Sigma}\) est une matrice de variance-covariance qui permet beaucoup plus de flexibilité en terme de variance et de dépendance entre les observations. En général, des paramètres seront associés à la matrice \(\mathbf{\Sigma}\) pour permettre différentes situations. Dans le cas d’un modèle linéaire standard, on a \(\mathbf{\Sigma} = \mathbf{I}\sigma^2\).


Ex.: variance différente pour les niveaux d’un facteur, variance augmentant avec certaines valeurs de \(x\), autocorrélation temporelle ou spatiale, etc.


Ces différentes situations peuvent être obtenues entre autres avec les packages nlme (fonctions gls ou lme) ou glmmTMB.

Modèles linéaires généralisés (GLM)

\[y \sim \mathcal{ExpFam}(\theta, ...)\]

\[\mathcal{g}(\mu)=\beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k\]

Structure d’un GLM

\[y \sim \mathcal{ExpFam}(\theta, ...)\]


\[ \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k\]


\[ g(E(y)) = g(\mu) = \eta \]

Composante systématique

\[ \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k\]

Le prédicteur linéaire \(\eta\) permet de relier les variables explicatives à la variable réponse.

Composante aléatoire

\[ y \sim \mathcal{ExpFam}(\theta) \]


Les observations \(y\) suivent une distribution particulière appartenant à la famille des distributions exponentielles avec des paramètres \(\theta\). La distributions normale, de Poisson, Binomiale et Gamma font parties de la famille des distributions exponentielles.

Fonction de lien

\[ g(E(y)) = g(\mu) = \eta \]


La fonction de lien (link function) est une transformation qui permet de relier la réponse attendue \(E(y)\) (i.e. la moyenne) au prédicteur linéaire.


C’est ce qui permet de relier les variables explicatives à la variable réponse. Pour chaque distribution qui pourrait être utilisée dans le cadre d’un GLM, une fonction de lien dite canonique est utilisée par défaut, mais il existe plusieurs fonctions de liens possibles.

GLM Poisson

\[ y \sim \mathcal{Pois}(\lambda) \]


\[ log(\lambda) = \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k \]


On a donc des observations \(y\) qui suivent une distribution de Poisson avec un paramètre \(\lambda\). Une fonction de lien \(log\) permet de relier linéairement les variables explicatives à la valeur attendue \(E(y) = \lambda\) de la variable réponse.

Distribution de Poisson

Cette distribution est souvent utilisée pour décrire le nombre de fois qu’un événement survient dans un intervalle de temps ou d’espace.

\[f(k,\lambda) = Pr(y = k) = \frac{\lambda^{k}e^{-\lambda}}{k!}\]

Cette fonction permet de calculer la probabilité d’observer un événement \(k\) fois avec une distribution de Poisson ayant un paramètre \(\lambda\). La moyenne de cette distribution est égale à \(\lambda\).

Cette distribution a la propriété que \(\lambda = E(y) = Var(y)\). Autrement dit, la moyenne et la variance sont égales et correspondent au paramètre \(\lambda\) de la distribution.

D’où provient-elle?

Simulons une distribution de Poisson en distribuant aléatoirement des observations dans un carré 10 x 10. Ensuite, divisons ce carré en 100 cellules.

library(sf)

n<-200
x<-runif(n,0,10)
y<-runif(n,0,10)

sfc <- st_sfc(st_polygon(list(rbind(c(0,0), c(10,0), c(10,10), c(0,0)))))
grid <- st_make_grid(sfc, cellsize = 1)
plot(grid,axes=TRUE)
points(x,y)

Nb d’observations par cellule

Ensuite, calculons le nombre d’observations par cellule et illustrons la distribution du nombre d’observations pour chaque cellule. En moyenne, on devrait retrouver 2 observations par cellule (200 obs. / 100 cell. = 2).

o<-st_intersects(grid,st_as_sf(data.frame(x,y),coords=c("x","y")))
co<-sapply(o,length)
hist(co,breaks=seq(-0.5,max(co)+0.5,by=1),freq=FALSE,main="",xlab="Nb d'observations")
lines(0:max(co),dpois(0:max(co),lambda=2),type="b",cex=2)
legend("topright",lty=c(NA,1),legend=c("Observée","Théorique"),bty="n",cex=1,pch=c(22,21),pt.bg=c("gray","white"))

Ce serait la même chose si on simulait des observations aléatoires dans le temps et si on comptait le nb d’observations par intervalles de temps.

Simulons un modèle

Simulons un modèle à partir d’une distribution de Poisson.

n<-100
x<-runif(n,0,10)
lambda<-exp(-2+0.4*x)
y<-rpois(n,lambda)
d<-data.frame(y,x)

m<-glm(y~x,family="poisson",data=d)
summary(m)
## 
## Call:
## glm(formula = y ~ x, family = "poisson", data = d)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.46287    0.30996  -7.946 1.93e-15 ***
## x            0.45294    0.03784  11.971  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 289.290  on 99  degrees of freedom
## Residual deviance:  85.894  on 98  degrees of freedom
## AIC: 245.83
## 
## Number of Fisher Scoring iterations: 5

Coefficients

Les coefficients estimés sont très similaires à ce qui a été utilisé pour générer le modèle (-2 et 0.4).

coef(m)
## (Intercept)           x 
##  -2.4628709   0.4529379

Prédictions

Illustrons l’effet de la variable explicative.

visreg(m)

Sous quelle échelle?

visreg(m,scale="response")

Échelle du prédicteur linéaire ou de la réponse

par(mfrow=c(1,2))
visreg(m)
visreg(m,scale="response")

Remarquez que dans le graphique de gauche, les valeurs en y sont négatives ou positives. En fait, la fonction log permet de transformer les comptes allant de 0 à l’\(+\infty\) en valeurs allant de \(-\infty\) à \(+\infty\). Sous l’échelle log, les valeurs entre 0 et 1 sont négatives et les valeurs supérieures à 1 sont positives.

\[ log(1) = 0 \] \[ log(0.1) = -2.302586 \] \[ log(5) = 1.609438 \]

ggeffects

g<-ggpredict(m,terms=c("x [n=50]"))
plot(g,add.data=TRUE,jitter=FALSE)

predict

v<-seq(min(x),max(x),length.out=50)
plot(y~x,data=d)
p<-predict(m,data.frame(x=v))
lines(v,p)

Argument type

Par défaut, l’argument type = "link" et il faut spécifier type = "response" pour obtenir les prédictions sous l’échelle de la réponse.

v<-seq(min(x),max(x),length.out=50)
plot(y~x,data=d)
p<-predict(m,data.frame(x=v),type="response")
lines(v,p)

Suppositions de base

On pourrait être tenté d’utiliser la fonction plot pour vérifier les conditions d’applications.

par(mfrow=c(2,2),mar=c(4,4,2,2))
plot(m)

Quelles sont les suppositions de bases?

\[ y \sim \mathcal{Pois}(\lambda) \]

\[ log(\lambda) = \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k \]

Certaines suppositions sont les mêmes que pour un modèle linéaire standard avec une distribution, mais d’autres diffèrent.

Vérification avec DHARMa

Une façon de vérifier les supposition de bases des modèles est d’utliser le modèle pour simuler des observations et ensuite comparer ces observations simulées aux observations réelles. Si les observations réelles ne ressemblent pas aux observations simulées à partir du modèle, c’est que le modèle ne représente pas bien ce qui se passe avec les données réelles et donc que qqch devrait être modifié dans le modèle.

Le package DHARMa se base sur ces simulations pour comparer les observations attendues par le modèle aux observations réelles. Plus spécifiquement, la méthode utilise les résidus quantiles pour comparer chaque observation aux observations attendues pour cette observations. Les résidus quantiles sont rapportés entre 0 et 1 et si le modèle est approprié pour les données, on devrait s’attendre à ce que les résidus soient distributés uniformément entre 0 et 1 pour toutes les observations.

library(DHARMa)

s<-simulateResiduals(m)
plot(s)

Si le modèle est appropré pour les données, on devrait voir ici une distribution uniforme des points et il ne devrait pas y avoir de patrons particuliers dans la distribution des points dans le graphique de droite. En particulier, les lignes de prédictions illustrées devraient être alignées sur les trois lignes représentant les quantiles 25, 50 et 75%. Dans le graphique quantile-quantile de gauche, les points devraient être alignés sur la droite.

Autres vérifications

Les problèmes de sous- ou sur- dispersion et de surabondance de 0 (zéro-inflation) sont très fréquents avec les GLM. Ces deux fonctions de DHARMa permettent d’évaluer si les observations réelles diffèrent des observations simulées en terme de présence de 0 ou de disperison. En rouge, ce sont les observations réelles et en gris les observations simulées. Si les observations réelles sont bien en dehors des observations simulées, C’est qu’il y a un problème avec le modèle et il faut un ajustment.

par(mfrow=c(1,2))
testDispersion(s)
testZeroInflation(s)

GLM binomial

\[ y \sim \mathcal{Binom}(n,p) \]

\[ logit(p) = log\left(\frac{p}{1-p}\right) = \eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_k x_k \]

Régression logistique

\[ y \sim \mathcal{Binom}(1,p) = \mathcal{Bernoulli}(p)\]

La distribution de Bernoulli est un cas particulier de la distribution binomiale avec un tirage (e.g. un individu est vivant ou mort, un individu est mâle ou femelle, etc.)

Distribution binomiale

\[f(k,n,p) = Pr(y = k) = \binom{n}{k}p^k(1-p)^{n-k}\]

\[\binom{n}{k}=\frac{n!}{k!(n-k)!}\]


Ce qui représente le nombre de façons de prendre \(k\) éléments parmi \(n\) éléments.


Cette fonction permet de calculer la probabilité d’observer \(k\) succès si \(n\) tirages sont faits et que la probabilité de succès est de \(p\)


Cette distribution a la propriété que \(E(y) = np\) et \(Var(y) = np(1-p)\) , autrement dit, la moyenne et la variance sont entièrement détemrinées par la propbabilité de succès et par le nombre de tirages. On a donc une autre distribution assez restrictive.

Exemple

Supposons qu’on tire pile (0) ou face (1) 10 fois et qu’on calcule le nombre de face et qu’on répète cette exercice 500 fois. Quelle sera la distribution du nombre de faces obtenus?

n<-500
co<-sapply(1:n,function(i){
  s<-sample(0:1,10,prob=c(0.5,0.5),replace=TRUE)
  sum(s) # on additione les faces (1)
})

hist(co,breaks=seq(-0.5,max(co)+0.5,by=1),freq=FALSE,main="",xlab="Nb de faces")
lines(0:max(co),dbinom(0:max(co),prob=0.5,size=10),type="b",cex=2)
legend("topright",lty=c(NA,1),legend=c("Observée","Théorique"),bty="n",cex=1,pch=c(22,21),pt.bg=c("gray","white"))

Simulons un modèle

logit<-function(p){
  log(p/(1-p))
}

inv.logit<-function(x){
   exp(x)/(1+exp(x))  
}

n<-200
x<-runif(n,0,10)
z<--2+0.5*x
y<-rbinom(n,size=1,prob=inv.logit(z))

d<-data.frame(y,x)

m<-glm(y~x,family="binomial",data=d)
summary(m)
## 
## Call:
## glm(formula = y ~ x, family = "binomial", data = d)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.68335    0.43063  -6.231 4.63e-10 ***
## x            0.63864    0.08853   7.214 5.45e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 272.74  on 199  degrees of freedom
## Residual deviance: 185.42  on 198  degrees of freedom
## AIC: 189.42
## 
## Number of Fisher Scoring iterations: 5

Visualisation

Avec visreg

par(mfrow=c(1,2))
visreg(m)
visreg(m,scale="response")

Remarquez que dans le graphique de gauche, les valeurs en y sont négatives ou positives. En fait, la fonction logit permet de transformer les probabilités allant de 0 à 1 en valeurs allant de \(-\infty\) à \(+\infty\). Sous l’échelle logit, une valeur de 0 corrrespond à une probabilité de 0.5.

\[ log\left(\frac{p}{1-p}\right) = log\left(\frac{0.5}{1-0.5}\right) = log(1) = 0 \]

Visualisation

Avec ggeffects

g<-ggpredict(m,terms=c("x [n=50]"))
plot(g,add.data=TRUE)

Suppositions de base

par(mfrow=c(2,2))
s<-simulateResiduals(m)
testUniformity(s)
testQuantiles(s)
testDispersion(s)
testZeroInflation(s)

Problèmes

Les problèmes les plus fréquents sont la surdispersion et la surabondance de zéros (zero-inflation)

Surdispersion

C’est lorsque les données sont plus variables que ce qui est attendu par le modèle. Les distribution de Poisson et binomiale sont très restrictives par rapport à leur variance et les problèmes de surdispersion sont très fréquents.


En général, une surdisperion ignorée fera en sorte que l’incertitude sera sous-estimée (trop de confiance dans nos résultats).


Plusieurs causes possibles (variables explicatives manquantes, phénomènes d’aggrégation, hétérogénéité, etc.)


Solutions:

Poisson -> Négative binomiale (glm.nb, glmmTMB, brms, etc.)

Binomiale -> Béta binomiale (glmmTMB, brms, etc.)

Obervation-Level Random Effects (OLRE) (voir Harrison 2015 et Harrison 2014).

Exemple

Voici un exemple de GLM de Poisson avec de la surdispersion. Pour créer cette surdispersion, j’utilise une distribution négative binomiale pour générer les observations au lieu d’une distribution de Poisson.

n<-100
x<-runif(n,0,10)
lambda<-exp(-2+0.4*x)
y<-rnbinom(n,mu=lambda,size=0.5) # 
d<-data.frame(y,x)

m<-glm(y~x,family="poisson",data=d)
g<-ggpredict(m,terms="x")
plot(g,add.data=TRUE)

Vérification avec DHARMa

Les résidus ne sont clairement pas distributés uniformément sur le second graphique et la surdispersion est évidente sur le 3e graphique.

s<-simulateResiduals(m)
par(mfrow=c(1,4))
testUniformity(s)
testQuantiles(s)
testDispersion(s)
testZeroInflation(s)

Correction avec glmmTMB

library(glmmTMB)

m<-glmmTMB(y~x,family="nbinom1",data=d)

s<-simulateResiduals(m)

par(mfrow=c(1,4))
testUniformity(s)
testQuantiles(s)
testDispersion(s)
testZeroInflation(s)

Sousdisperion

C’est lorsque les données sont moins variables que ce qui est attendu par le modèle. Ceci est beaucoup plus rare et probablement moins dommageable (e.g. taille de couvée/portée chez les oiseaux/mammifères voir Brooks et al. 2019).


En général, une sous-dispersion ignorée fera en sorte que l’incertitude sera surestimée (trop conservateurs dans nos résultats) ce qui est un moindre mal.


Solutions:

Poisson -> glmmTMB, family = "compois"

Zéro-inflation

Surabondance de 0 dans les observations par rapport à ce qui est attendu par le modèle. Il est possible aussi d’avoir une sous-abondance de 0 et/ou encore une surabondance de 1, mais ce dernier cas est un peu plus rare.


Solutions:

modèle zero-inflated (glmmTMB,pscl,brms, etc.)


À noter qu’une surdipersion peut être due à une surabondance de 0 et vice-versa. Ce sont deux problèmes qui peuvent interagir.

Extensions


Gamma    \(]0,\infty\)

Variable réponse continue, strictement positive et avec une asymétrie à droite.

Souvent une alternative à une transformation \(log\).

Ex.: concentration d’un composé, biomasse, etc.


Tweedie    \([0,\infty\)

Variable réponse continue positive et incluant des 0.

Ex.: biomasse capturée, précipitations, etc.

Voir Lecomte et al. 2013


Beta    \(]0,1[\)   ou   \([0,1]\)

Variable continue entre 0 et 1 exclusivement.

Pour proportion ou % ne provenant pas d’un ratio entre entiers.

Si 0 et 1 possibles, certaines transformations peuvent être utilisées.

Ex.: proportion d’une surface affectée, proportion de temps passé en alimentation, etc.

Voir Douma & Weedon 2019

2. Modèles mixtes/hiérarchiques



\[ y = \alpha + \beta x + \epsilon \]

\[ \epsilon \sim \text{Normale}(0,\sigma^{2}) \]

ou


\[ \mu = \alpha + \beta x \]

\[ y \sim \text{Normale}(\mu,\sigma^{2}) \]

Modèle linéaire hiérarchique simple



\[ y_{ij} = \alpha + \beta x_{ij} + \tau_{i} + \epsilon_{ij}\]

\[ \tau_{i} \sim \text{Normale}(0,\,\sigma_{\tau}^{2}) \]

\[ \epsilon_{ij} \sim \text{Normale}(0,\sigma_{\epsilon}^{2})\ \]



où il y a \(i\) groupes et \(j\) observations par groupe

Autre formulation



\[ \mu_{ij} = \alpha_{i} + \beta x_{ij}\]

\[ \alpha_{i} \sim \text{Normale}(\alpha,\,\sigma_{\alpha}^{2})\ \]

\[ y_{ij} \sim \text{Normale}(\mu_{ij},\,\sigma_{\epsilon}^{2})\ \]

Un modèle est dit hiérarchique lorsque certains de ses paramètres sont définis par une probabilité de distribution.

Simulons un exemple

set.seed(12345)
ng <- 20
nobs <- 10
n <- ng * nobs
id <- gl(ng, nobs)

# valeurs fictives des variables explicatives selon une distribution uniforme
x <- runif(n, 0, 100)

# valeurs des paramètres
a <- 100
b <- 1

# on génère l'effet aléatoire pour chaque groupe
rea <- rep(rnorm(ng, 0, 20), each = nobs)

# prédicteur linéaire 
mu <- (a + rea) + b * x

# observations à partir d'une distribution normale avec un écart-type de 5
y <- rnorm(n, mu, 5)

# données
d <- data.frame(y, x)

Illustrons le modèle

library(lme4)
library(ggeffects)

m <- lmer(y ~ x + (1|id), data = d)

g <- ggpredict(m, terms = c("x", "id [sample=20]"), type = "random")
plot(g, add.data = TRUE, show_ci = FALSE, colors = rep("black", 20), show_legend = FALSE)

Illustrons les pentes par individus

g <- ggpredict(m, terms = c("x", "id [sample=20]"), type = "random")
plot(g, add.data = TRUE, show_ci = FALSE, colors = sample(colors(), 20))

Retrouvons les paramètres

La variance de l’effet aléatoire était de 20² et la variance résiduelle était de 5².

summary(m)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ x + (1 | id)
##    Data: d
## 
## REML criterion at convergence: 1314.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.50905 -0.67270 -0.04039  0.61012  2.58655 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  id       (Intercept) 386.9    19.67   
##  Residual              25.3     5.03   
## Number of obs: 200, groups:  id, 20
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 101.35098    4.46512   22.70
## x             0.99676    0.01262   78.97
## 
## Correlation of Fixed Effects:
##   (Intr)
## x -0.153

Récupérons ces effets

fixef(m)
## (Intercept)           x 
## 101.3509752   0.9967622
ranef(m)
## $id
##    (Intercept)
## 1    2.8823702
## 2  -22.7874952
## 3    5.4199340
## 4  -28.1373544
## 5    1.0787346
## 6  -10.6043879
## 7   -6.4207965
## 8   31.3275567
## 9  -10.5628148
## 10   4.4707573
## 11 -24.5843971
## 12 -26.8974001
## 13  24.9063468
## 14  23.4502635
## 15  -2.0989706
## 16 -12.2974161
## 17  -1.4157682
## 18   8.9803887
## 19  44.1496752
## 20  -0.8592262
## 
## with conditional variances for "id"
apply(ranef(m)$id, 2, sd)
## (Intercept) 
##    19.60505

Coefficients

coef(m)
## $id
##    (Intercept)         x
## 1    104.23335 0.9967622
## 2     78.56348 0.9967622
## 3    106.77091 0.9967622
## 4     73.21362 0.9967622
## 5    102.42971 0.9967622
## 6     90.74659 0.9967622
## 7     94.93018 0.9967622
## 8    132.67853 0.9967622
## 9     90.78816 0.9967622
## 10   105.82173 0.9967622
## 11    76.76658 0.9967622
## 12    74.45358 0.9967622
## 13   126.25732 0.9967622
## 14   124.80124 0.9967622
## 15    99.25200 0.9967622
## 16    89.05356 0.9967622
## 17    99.93521 0.9967622
## 18   110.33136 0.9967622
## 19   145.50065 0.9967622
## 20   100.49175 0.9967622
## 
## attr(,"class")
## [1] "coef.mer"

Complexifions le modèle

b
## [1] 1
reb <- rep(rnorm(ng, 0, 0.5), each = nobs)

mu <- (a + rea) + (b + reb) * x

y <- rnorm(n, mu, 5)

d <- data.frame(y, x)

Illustrons les pentes

m <- lmer(y ~ x + (x|id), data = d)

g <- ggpredict(m, terms = c("x", "id [sample=20]"), type = "random")
plot(g, add.data = TRUE, show_ci = FALSE, colors = rep("black",20), show_legend = FALSE)

Retrouvons les paramètres

La variance de l’effet aléatoire sur la pente était de (0.5)².

summary(m)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ x + (x | id)
##    Data: d
## 
## REML criterion at convergence: 1370.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.19310 -0.58325  0.02874  0.58522  2.42595 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  id       (Intercept) 388.2872 19.7050       
##           x             0.2105  0.4588  -0.27
##  Residual              22.0018  4.6906       
## Number of obs: 200, groups:  id, 20
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 100.4564     4.4699  22.474
## x             0.8959     0.1033   8.672
## 
## Correlation of Fixed Effects:
##   (Intr)
## x -0.283
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00784532 (tol = 0.002, component 1)

Récupérons ces effets

fixef(m)
## (Intercept)           x 
## 100.4563555   0.8958772
ranef(m)
## $id
##    (Intercept)           x
## 1    -4.259743 -0.11207853
## 2   -24.555894  0.09388618
## 3     9.760813  0.08930651
## 4   -23.972706 -0.43249104
## 5     6.379230  0.06853379
## 6   -10.672292  0.92265805
## 7    -4.656865 -0.14388477
## 8    33.097624 -0.22946144
## 9    -7.593899  0.23044679
## 10   -2.310431  0.38486838
## 11  -32.815439 -0.47826203
## 12  -26.583571  0.53612804
## 13   22.768061  0.34741573
## 14   22.280452 -0.58389607
## 15    1.432798  0.75717249
## 16   -8.201353 -0.31288440
## 17    1.063482  0.10990165
## 18   10.101182  0.03586551
## 19   39.745709 -0.88012487
## 20   -1.007156 -0.40309995
## 
## with conditional variances for "id"
apply(ranef(m)$id, 2, sd)
## (Intercept)           x 
##  19.4359900   0.4557504
coef(m)
## $id
##    (Intercept)          x
## 1     96.19661 0.78379868
## 2     75.90046 0.98976339
## 3    110.21717 0.98518372
## 4     76.48365 0.46338617
## 5    106.83559 0.96441100
## 6     89.78406 1.81853525
## 7     95.79949 0.75199243
## 8    133.55398 0.66641577
## 9     92.86246 1.12632400
## 10    98.14592 1.28074559
## 11    67.64092 0.41761518
## 12    73.87278 1.43200525
## 13   123.22442 1.24329294
## 14   122.73681 0.31198114
## 15   101.88915 1.65304970
## 16    92.25500 0.58299281
## 17   101.51984 1.00577886
## 18   110.55754 0.93174272
## 19   140.20206 0.01575234
## 20    99.44920 0.49277726
## 
## attr(,"class")
## [1] "coef.mer"

Formulations

m <- lmer(y ~ x + (1|id), data = d)

m <- lmer(y ~ x + (x|id), data = d) # équivalent à (1 + x|id)

m <- lmer(y ~ x + (0 + x|id), data = d)

m <- lmer(y ~ x + (1|id) + (0 + x|id), data = d)

Corrélations intercepts/pentes

m <- lmer(y ~ x + (1|id) + (0 + x|id), data = d)


Le dernier est équivalent à

m <- lmer(y ~ x + (x||id), data = d)


et signifie qu’on assume que l’intercept et la pente ne sont pas corrélées ou qu’on ne veut pas estimer cette corrélation. En pratique, il n’y a pas énormément de raison de faire cela et par défaut le modèle va estimer cette corrélation.

Effets aléatoires croisés

Supposons que nous avons des individus observés plusieurs fois par année et ce sur plusieurs années.

m <- lmer(y ~ x + (1|id) + (1|year), data = d)

Effets aléatoires nichés

Supposons maintenant que les individus sont annuels et qu’ils ne sont observés que pour des années particulières.

m <- lmer(y ~ x + (1|year/id), data = d)


Équivalent à:

m <- lmer(y ~ x + (1|year) + (1|year:id), data = d)

Quand utiliser les effets aléatoires?


  • Mesures répétées sur certaines unités ou individus

  • Les observations sont regroupées sous des structures ou ne sont pas indépendantes

  • Il y a de la réplication à plusieurs niveaux de la hiérarchie

  • On veut estimer la variance ou généraliser à une population d’unités d’échantillonnage

Formulation matricielle



\[\boldsymbol{y} \sim \text{Normale}(\mathbf{\beta}\mathbf{X}, \mathbf{I}\sigma^2)\]

\[\boldsymbol{y} \sim \text{Normale}(\mathbf{\beta}\mathbf{X}, \mathbf{\Sigma})\]

Effet spatial, temporel, phylogénétique


Autocorrélation spatiale, temporelle, phylogénétique, etc…

Régularisation (shrinkage)


Complete pooling

No pooling

Partial pooling

Problèmes, questionnements



La GLMM FAQ!


https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html